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The river bar instability is revisited, using a liydrodynamical model based on Reynolds 
averaged Navier-Stokes equations. The results are contrasted with the standard analysis 
■ based on shallow water Saint- Venant equations. We first show that the stability of both 

transverse modes (ripples) and of small wavelength inclined modes (bars) predicted by 
'"O ■ the Saint- Venant approach are artefacts of this hydrodynamical approximation. When 

using a more reliable hydrodynamical model, the dispersion relation does not present 
any maximum of the growth rate when the sediment transport is assumed to be locally 
saturated. The analysis therefore reveals the fundamental importance of the relaxation 
of sediment transport towards equilibrium as it it is responsible for the stabilisation of 
small wavelength modes. This dynamical mechanism is characterised by the saturation 
number, defined as the ratio of the saturation length to the water depth Lsat/H. This 
, dimensionless number controls the transition from ripples (transverse patterns) at small 

Q-f Lsa.t/H to bars (inclined patterns) at large LsaxjU. At a given value of the saturation 

number, the instability presents a threshold and a convective-absolute transition, both 
controlled by the channel aspect ratio /3. We have investigated the characteristics of 
the most unstable mode as a function of the main parameters, Lsa,t/H, /? and of a 
subdominant parameter controlling the relative influence of drag and gravity on sediment 
transport. As previously found, the transition from alternate bars to multiple bars is 
mostly controlled by the river aspect ratio /?. By contrast, in the alternate bar regime 
(large Lsat/ H), the selected wavelength does not depend much on /3 and approximately 
scales as H'^/'^L\1^ /C , where C is the Chezy number. 



1. Introduction 

Bars are large scale river bedforms whose wavelength is comparable to channel width 
and whose height is comparable to flow depth. Due to these two large length-scales, bars 
interact with structures and navigation, and have a profound impact on several aspects of 
river engineering like river control (channel shift) and hazard prevention (bairk failures) . 

Free bars form in straight or weakly curved channels and consist of repetitive sequences 
of migrating pools and diagonal depositional riffles. Alternate (single row) bars are the 
generic type of free bars in sandy streams and gravel bed rivers, when the channel is 
narrow enough (Callander, 1969). When the channel presents a large width to depth 
aspect ratio (typically in glaciary valleys), the river develops multiple bars, leading to 
braiding patterns (Fujita & Muramoto, 1985; Crosato & Mosselman, 2009). Fluvial bars 
may also be forced by various effects like curvature, width variations or confluences. The 
formation of free alternate bars in an initially straight channel has been considered as 
the possible origin of incipient meandering (Fukuoka, 1989). The 'bar theory' of river 
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meandering implies the coexistence of free and forced bars in weakly curved channels 
and transition from migrating free to steady forced bars in developing meanders (Fujita 
and Muramoto, 1982; Tubino & Scminara, 1990; Crosato et al., 2011). This idea has 
been further developed theoretically in the framework of the so-called 'bend theory', in 
which the formation of meanders results from the non-linear interaction between self- 
excited (free) and curvature-driven (forced) bed responses (Seminara & Tubino, 1989). 
Overall, the emergence of bars and their coupling with the shape of the banks are crucial 
ingredients required to understand the global morphodynamics of alluvial channels at 
the river scale (see the recent review of Seminara, 2010). 

Since the seminal work of Kennedy (1963), there is a wide agreement that bedforms 
- both free bars and sand ripples - develop by linear instability of an erodible bed subject 
to a turbulent flow: infinitesimal perturbations of the bed (and therefore of the flow) are 
selectively amplified at certain wavelengths in the course of their downstream migration. 
By contrast with sand ripples, which are transverse to the fiow, bars develop from a 
three-dimensional instability of the sediment bed in which the most unstable modes are 
inclined with respect to the flow (Chang & Simons, 1970; Chang ct al., 1971; Schumm 
& Khan, 1972; Ikeda, 1983; Fujita & Muramoto, 1985; Lisle et al., 1991, 1997; Lanzoni, 
2000a, 2000b; Devauchelle et al., 2010a, 2010b; Andreotti et al., 2012). Like optical 
or acoustical wave-guides do, the channel selects a discrete number of guided unstable 
modes labelled by the number m of rows of bar perturbations (to = 1 corresponds to 
alternate bars). The transverse wave- number is therefore directly selected by the chan- 
nel width (Parker, 1976). Several increasingly refined linear stability analysis have been 
published in the literature (Callander, 1969; Engelund & Skovgaard, 1973; Parker, 1976; 
Freds0e 1978; Colombini et al., 1987; Garcia & Niiio, 1993). These theories predict the 
dispersion relation - both the growth rate a^k^) and the propagation velocity c{kx) of 
perturbations of wavcnumbcr - within the linear regime, the marginal stability condi- 
tions in the space of flow and sediment parameters and the wavenumber fcmax selected by 
the instability mechanism, i.e. those corresponding to the maximum growth rate. The 
parameter controlling the instability is the width to depth ratio /? = W/H characterising 
the cross section of the channel. Colombini et al. (1987) have shown that nonlinear 
effects cause bed perturbations to reach an equilibrium flnite amplitude characterized 
by periodic diagonal fronts. Under suitable conditions, which are rarely encountered in 
nature, the above equilibrium solution may in turn become unstable and bifurcate into a 
more complex quasi-periodic pattern (Schielen et al., 1993). One of the most important 
property of the linear instability is the existence of a particular value of /3 for which 
the migration speed of the marginally stable mode vanishes (Blondeaux & Seminara, 
1985). This so-called 'resonant' condition is the key ingredient of the 'bend theory' of 
meandering, which assumes a weak coupling between bars and banks. River meanders 
are assumed to excite the bar instability without affecting the shape of the modes nor 
the dispersion relation; then, the coupling is optimal when non-growing free bars do not 
propagate with respect to the banks. 

Recent reviews have claimed that the topic can be considered as fairly settled (Tubino 
et al., 1999; Seminara, 2010). However, most existing theories (all but Engelund & 
Skovgaard (1973)) are based on depth-averaged Saint- Venant shallow water equations, 
which assume that the length-scale over which the flow varies is much smaller than 
the flow thickness H. Our primary goal, here, is to investigate the robustness of the 
results previously obtained, when the hydrodynamical disturbances are computed with 
Reynolds averaged Navier Stokes equations, which describe the vertical structure of the 
flow. Moreover, most of the previous papers do not take into account the fact that 
sediment transport adapts to a variation of the flow strength with a spatial lag. This 
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sediment transport saturation transient, characterised by the so-called saturation length 
Lsat , plays a dominant role in the linear instability of subaqueous sand ripples and aeolian 
dunes, as it controls the emergent wavelength (Andreotti ct al., 2002; Lagree, 2003; 
Elbclrhiti et al., 2005; Valance, 2005; Valance & Langlois, 2005; Claudin & Andreotti, 
2006; Charru, 2006; Fourriere et al., 2010; Duran et al., 2011). The saturation length 
allows one to build a new dimensionless number, the saturation number Ls^t/H, which 
compares the stabilising effect of the sediment transport relaxation transient to that of 
the free surface (Andreotti et al., 2012). Our second aim is to look at the influence of 
this dimensionless parameter on the bar instability. 

In section[5J we focus on three dimensional hydrodynamics above a modulated bottom, 
comparing the results obtained with Saint- Venant shallow water (SVSW) equations and 
with Reynolds averaged Navier-Stokes (RANS) equations. We show that shallow water 
approximation agrees semi-quantitatively with the full hydrodynamical resolution pro- 
vided that the wavelengths in streamwise and perpendicular directions are larger than 
the flow thickness rescaled by the Chezy number. At small wavelength, however, they are 
quantitatively (wrong orders of magnitude) and qualitatively (wrong signs) incorrect. In 
section[31 we investigate the consequences of these discrepancies on the bar linear stability 
analysis. We show that the maximum growth rate found in previous papers is an artefact 
of Saint- Venant equations which disappears when using three-dimensional hydrodynam- 
ics: large wavenumbers are incorrectly predicted to be stable. As a consequence, the 
introduction of the saturation length igat is necessary to recover an instability forming 
bars. We then show that the emergence of bedforms in a channel is controlled primarily 
by two parameters: the width to depth ratio /? = W/H, as found in former studies, but 
also the saturation length rescaled by the flow thickness Lsat/H. In section |4l we focus 
on the RANS based approach to revisit the basic properties of the bar instability, and 
investigate in particular the influence of the saturation parameter Lsat/H as well as the 
fluid shear velocity. We then study the convective-absolute transition of the instability. 
We finally present the stability diagram showing the transitions from ripples to alternate 
bars and from from alternate bars to multiple bars. 

2. Hydrodynamics over bedforms: contrasting Saint- Venant 
approximation with a Reynolds averaged calculation 

In this section we focus on the purely hydrodynamical part of the problem. We in- 
troduce Saint- Venant shallow water and Reynolds averaged Navier-Stokes equations and 
compare their predictions for the modulation of the basal shear stress over a gently 
undulated bottom. 



We consider a turbulent stream whose width W is larger than the flow depth H. It 
flows on a plane inclined at an angle 6 with respect to the horizontal. As schematised in 
Fig- HI X is the direction of the flow, y is in-plane transverse and z is normal to the plane, 
oriented upwards. We do not describe the lateral boundary layers, of thickness ~ H and 
assume that the banks act as frictionless impenetrable boundaries. Then, the base flow is 
similar to the homogeneous situation - unbounded in the x and y directions ~ for which 
the pressure profile p{z) is hydrostatic and the velocity profile Ux{z) is well described by 
a logarithmic law (Julien, 1998): 



2.1. Geometry and base flow 




(2.1) 
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Figure 1. (a) Free alternate bars in the Loire river, close to St-Mathurin-sur-Loire (France; 
47° 24' 50" N, 0° 19' 04" W). Photo credit: S. Rodrigues, universite Frangois Rabelais of Tours. 
These bars, of wavelength A ~ 2000 m, mostly evolve during annual floods, characterised by a 
typical flow depth // ~ 4 m. Sediment transport then mostly takes place in suspension, with 
a saturation length Lsat around 160 m. The saturation number Lsax/H is therefore around 40. 
As the river width W is around 480 m, the aspect ratio j3 — W/H is around 120. (b) Multiple 
bars in the Congo river, with a similar mean wavelength around A ~ 2000 m. The relevant flow 
depth H is uneasy to define but is smaller than 10 m. As the river width W is around 9000 m, 
the typical aspect ratio /? = W/H lies in the range 1000 — 2000. (c) Schematic showing a plane 
wave mode inclined by an angle a with respect to the flow direction, (d) Schematic showing 
free alternate bars in a channel of width W. They correspond to guided modes obtained by 
superimposing two plane waves propagating with opposite angles. The notations are the same 
as in Andreotti et al. (2012). 
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z = is the location of the bottom, where the velocity vanishes, g is the gravity 
acceleration; p is the fluid density, assumed constant; = \/gH sin 9 is the basal shear 
velocity; k = 0.4 is the von Karman phenomcnological constant; zq is the hydrodynamical 
bottom roughness and is in general much smaller than H. Throughout this paper, we 
will use H/zq = 10^, which is a typical value for natural rivers. The depth averaged 
velocity u = jj Ux{z)dz is related to the shear velocity by the Chezy coefficient C: 



C 



where we have used the hypothesis that H/zq 3> 1. The y- and z-components of the 
velocity are null. We define the Froudc number F as the ratio of the surface velocity to 
the velocity of gravity surface waves in the shallow water approximation: 

T= ~ -Vsing In— . 2.4) 

\J gU K Zq 

This number can be of order unity in fiumcs but is in general small for large natural 
rivers, due to their small slopes. 

2.2. Linear response of a flow to steady bedforms 

Before addressing the formation of bars, we discuss the linear response of hydrodynam- 
ics to bed perturbations. The technical aspects of the present analysis are detailed in 
Fourricre et al. (2010) and Andreotti et al. (2012). Therefore, we shall here only sketch 
the main lines of the computations and refer the interested reader to these papers. As 
the base state of the bed is homogeneous in x and y, the bcdform elevation Z{x,y) can 
be decomposed over normal Fourier modes. Without loss of generality, we can therefore 
consider one such mode 

Z = - tan 61 a; + (e'''^^°' " ^+"'" " = - tan 61 a; + Ce'^'^^ (2.5) 

corresponding to undulations making an angle a with the direction of the flow (see 
Fig. HI), k = {kx,ky) = (fc cosa, fc sina) is the wave vector. We shall compute the first 
order linear response of the flow to this bed perturbation, describing hydrodynamics 
either by (i) depth averaged Saint- Venant equations, or (ii) Reynolds averaged Navier- 
Stokes equations. 

In the SVSW approach, the governing equations read in the limit of small slopes 

V • (hu) = 0, (2.6) 

{u-V)u^-gV{Z + h)-C^-^, (2.7) 

where u = {ux,Uy) is the two-dimensional depth averaged velocity and h is the local 
water depth. In the homogeneous case, these equations have the simple solution h = H 
and u = y^gH ta.n9/C ex ■ In this approximation, the basal shear stress f* is assumed 
to be a function of the velocity: = —Cp\u\u. Importantly, the typical length-scale in 
the plane (x, y) is H/C, which is much larger than the fiow depth H as a typical value 
of the Chezy number is C ~ 10^^. 

The RANS equations have exactly the same form as Navier-Stokes equations. 



djUj = 0, 



(2.8) 
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dtUi + UjdjUi = gi - -dip - -d^nj. (2.9) 

except that Tij is the Reynolds stress tensor, which includes both viscous and turbulent 
diffusion of momentum. Assuming a fully turbulent state, at large Reynolds number, 
we use here a first order Prandtl-like turbulent closure to relate to the strain rate 

^ pK^L^\i\ Qx'l7l % - 7.,) • (2.10) 

X is another phenomenological constant in the range 2.5-3. The mixing length is chosen 
equal to L = (z + zo)\/l — z/H, in order to recover the logarithmic flow profile for the 
base state fEg. 12. ip . By comparison to SVSW equations, there is one more parameter 
in RANS, the ratio H/ zq associated to the vertical structure of the flow. 

The bed undulation amplitude ^ is assumed to be vanishingly small in front of the 
flow depth and the wavelength, so that all equations can be linearised with respect to 
the small parameter kC,. First order disturbances are computed for the pressure, for the 
water depth and for all the components of the velocity and of the stress tensor. As shown 
below, in the context of sediment transport and bedform evolution, the central quantity 
is the basal shear stress t^^ (see section|3]). Denoting by a superscript A the space Fourier 
transform, the relationship of the basal stress linear response to the elevation proflle Z 
can be written as 

= -P"* (^x + iB..)kZ, and f^^, = -pul {Ay + iBy)kZ. (2.11) 

Ax, Bx, Ay and By are respectively the components of the basal shear stress in phase {A^ 
and Ay) and in quadrature (B^ and By) with the bed deformation. These coefficients 
have analytical expressions in the case of the Saint- Venant shallow water approximation, 
but are the result of numerical integration for the RANS equations. They are functions 
of kxH and kyH, and depend on T and H/zq. As developed in the next sections, 
the instability is controlled by the coefficients Bx and By and the bedform propagation 
speed by Ax and Ay. One shall therefore pay attention to their sign and to their order 
of magnitude. 

2.3. Basal shear stress coefficients 

As most former models of bar formation are based on SVSW equations, it is important to 
compare the behaviour of the basal shear stress coeffcients Ax, Bx, Ay and By obtained 
with this approximation, to the more reliable predictions of RANS. The stress coefficients 
are displayed in Figs. [2] and [3] as a function of kxH for a fixed transverse wavenumber 
kyH . The corresponding angle a = arctan ^ is therefore continuously decreasing from 
7r/2 in the limit /c^ji? — > to in the limit kxH — oo. 

Let us first focus on the case of a small transverse wavenumber kyH = 0.1 (Fig. [2]). 
SVSW and RANS predictions are qualitatively similar for in-pliase shear stress compo- 
nents: Ax starts negative at small kxH, changes sign and reaches a maximum around 
kxH = 0.1, while Ay is always positive, also with a maximum around the same wavenum- 
ber. The scaling laws obtained with SVSW and RANS are however different at large kxH. 
The structure of the flow disturbance above bedforms does not resemble the base profile 
but presents a layered structure originally discussed in Jackson & Hunt (1975). Such 
a structure cannot be described by the depth-averaged equations of the Saint- Venant 
approach. Moreover, the predictions for the components in quadrature are qualitatively 
different: RANS equations give an always positive Bx, while Saint- Venant equations 
predict that Bx changes sign and become negative at large kxH. RANS predicts, as 
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Figure 2. Discrepancies between hydrodynamical models at large k^H. (a-d) Saint-Venant 
shallow water (SVSW) and Reynolds averaged Navier-Stokes (RANS) predictions of the basal 
shear stresses coefficients Aj;, B^, Ay and By. They have been computed for a fixed transverse 
wavenumber kyH = 0.1 and are plotted against k^H. The other parameters are T = 0.1 and 
H/zq — 10^. Graphical convention for log-log plots: negative values are plotted in absolute 
value, but are displayed with a type of line different from positive ones (see legend), (e) Angle 
Q between k^ and ky as a function of k^H. In all panels, vertical grey dashed lines indicate 
the crossover from rather longitudinal patterns at small wavenumber k^H {a > n /A) to rather 
transverse patterns at large k^H {a < 7r/4). 



observed experimentally (see the review Charru et al., 2013), that the basal shear stress 
is phase advanced with respect to topography while SVSW predicts a phase delay. As 
this phase between flow and relief is precisely at the origin of the emergence of bcdforms, 
the consequences of this discrepancy are extremely important: the Saint-Venant approx- 
imation misses the instability at large wavenumber. For By, the situation is reversed: By 
computed from the RANS equations switches from negative values at small wavenumber 
to positive ones, again with a maximum around k^H = 0.1, while it stays negative for all 
kxH in Saint-Venant modelling. This disagreement becomes worse for larger transverse 
wavenumber kyH — 1 (Fig. |3]), with an additional quantitative difference: the coefficient 
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Figure 3. Discrepancies between hydrodynamical models at large kyH. Same as Fig. [J] but 

for kyH = 1. 



Ax now changes sign at wavcnumbers that differ by almost one decade from SVSW to 
RANS frameworks. 

These differences are not surprising: contrarily to the lubrication approximation in the 
viscous regime, Saint- Venant equations are not well-controlled approximation of Navier- 
Stokes equations. As they rely on the shallow water approximation, one expects them to 
provide predictions semi-quantitatively valid in the limit where k^H/C and kyH /C are 
both small. In the left sides of all panels of Fig. [2l when kxH is smaller than the Chezy 
number C ~ 10~^, the two hydrodynamical descriptions effectively match each other. 
However, Saint- Venant approximation leads to spurious results, including wrong signs for 
the most important shear stress coefficients, as soon as k^H/C or kyH/C are not small. 
Because we work here at fixed ky (see below), the SVSW equations become valid in the 
limit of inclined bedforms (large a), when kyH is smaller than C. We emphasise that 
the prediction of RANS equations for the basal shear stress disturbance is very robust to 
the choice of the turbulent closure: although the flow in the outer layer can be sensitive 
to this choice, the structure of the flow in the inner layer remains almost unchanged 
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(Ayotte et al., 1994; Charru et al., 2013). Although expernTiental confirmations of the 
predictions are sparse (Poggi et al., 2007; Claudin et al., 2012), the RANS predictions 
can be considered as reliable in the turbulent regime. 



3. Consequences of hydrodynamical modelling for the bar instability 

We now wish to investigate the consequences of the differences between SVSW and 
RANS hydrodynamical predictions on the bar instability. In this section we therefore per- 
form the linear stability analysis of a channelised flat erodible bed towards the formation 
of alternate bars. 



3.1. Sediment transport 

We start from a generic description of sediment transport, characterised by a volumetric 
flux q = {q^,q^) defined as the volume of grains - packed at the bed volume fraction 
- passing per unit time trough a vertical surface of unit width and infinite vertical 
extension. The value gsat of this flux in the homogeneous steady state is said 'saturated'. 
Both for bed load and suspended load, qsat is a function of the basal shear stress f* = 
i'^xz^'^yz) wefl as of the longitudinal and transverse slopes of the bed. We write it 
under the generic following form: 



Tth) t, 



(3.1) 



where is a dimensionfull constant of proportionality, Tth is the threshold value below 
which transport vanishes, and 7 is a semi-empirical exponent: 7 ~ 3/2 for bed-load and 
turbulent suspension, see e.g. Meyer-Peter & Miiller (1948), Einstein (1950), Fernandez 
Luque & van Beek (1976), Lajeunesse et al. (2010), Andreotti et al. (2012). The 
longitudinal slope of the bed modifies at the first order the threshold as 



Tth 



'th 



1 



-dxZ 



(3.2) 



where /i ~ 0.5 is the avalanche slope (Dey, 2003; Fernandez Luque & van Beek, 1976). 
Below we make use of the shear velocity threshold defined as u'^^ = yfr^Jp- The trans- 
verse slope comes into the expression of the unit vector i, which, at the first order, 
writes 



t 



(3.3) 



where e|| = /t^ is the unit vector in the direction of the basal shear stress and e±_ = 
{—Tyz^'''xz)l'^^ is its orthogonal counterpart (Andreotti et al. 2012). Combining the 
above expressions, the first order corrections to the saturation flux can be computed and 
read: 



Q 



Ax + iB^ - - [ — ] cos a 



kZ = Q{ax + ibx)kZ, 



ysat 



7 



^th 



-\- iBy 



-*th 



sma 



where we have deflncd the reference flux Q = 
flux response coefficients ax, bx, ay and by. 



717 



kZ = Q{ay + iby)kZ, 



(3.4) 
(3.5) 



N21T-1 



ui and introduced the 
Note that transverse transport component 



vanishes at the sediment transport threshold, as m* 



^th- 
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3.2. Dispersion relation, assuming a locally saturated sediment transport 

In the case of an erodible bed, we assume a time scale separation between the dynamics 
of the flow and that of the bed, which is more slower. The flow can thus be consid- 
ered as steady during the evolution of the bed profile, governed by the sediment mass 
conservation equation: 

dtZ + V-q^O. (3.6) 

We first assume, as usually done in former articles on the formation of bars, that sediment 
transport is in local equilibrium with the flow, i.e. q = ^sat evaluated with the local basal 
shear stress. We will see that this strong hypothesis must be relaxed and that it is 
necessary to take the relaxation of transport towards saturation into account. 

This linear stability analysis is performed to study bar formation, i.e. in the geometry 
of a channelized flow, with straight banks located at y = and y = W. We assume that 
the base profiles (|2.1[ 12. 2p are valid in most of the channel. As already stated, we expect 
the influence of the banks to be limited to a boundary layer whose thickness is on the 
order of H. We do not describe this region close to the banks in details, and impose a slip 
boundary condition: the velocity of the fluid at the bank must be tangent to the bank 
itself. To obtain a mode of instability guided by the channel, one superimposes plane 
waves of transverse wavenumbers ky and —ky. The condition of null transverse velocity 
at the banks selects a discrete number of transverse wavenumbers labelled by the mode 
number m: kyW = ttitt. Equivalently, this condition can be expressed in terms of the 
aspect ratio /? = W/H as: 

kyH = m ^ . (3.7) 

m = corresponds to transverse bedforms (sand ripples) and m = 1 corresponds to 
alternate bars. We will first discuss the modes inclined with respect to the flow for 
a fixed order m > 0. Then, we will summarise the results obtained for m = by 
Fourriere et al. 2010. We will finally compare the growth rate of the different orders 
m for the phase diagram of Fig. [T^ Defining the growth rate a and the phase velocity 
c in Z{x, y,t) = ( cos{kyy) exp[at + ik^ {x — ct)] , one obtains, from Eq. 13.61 and local 
equilibrium condition q = cfsat : the following dispersion relation: 

a — ik^c = ^ik Q [fc^, (a^ + ih^.) + ky [ay + ihy)] . (3-8) 

As can be seen from Eqs. 13.4 1 and 13.51 where the coefficients a^, h^, Oy and hy arc defined, 
both a and c depend on the ratio u^ju^^ - besides k^H, kyH, J- and H/zq. 

The growth rate cr and the phase velocity c are plotted in Fig. S] as functions of k^H, 
for a given set of parameters. One observes that the discrepancies between hydrody- 
namical predictions obtained with SVSW and RANS have fundamental consequences on 
the instability. A first difference affects the growth rate curves: while SVSW descrip- 
tion shows an unstable region (cr > 0) with a most unstable mode corresponding to an 
intermediate wavenumber fcmax^^^ — 0.05 (triangle in Fig. the RANS liydrodynam- 
ical calculation does not present any maximum growth rate for kxH of order one or 
smaller. This fundamental discrepancy results from the fact that Saint- Venant equations 
predict the wrong sign for for this range of k^H, i.e. a phase delay of the basal shear 
stress with respect to topography instead of a phase advance. Large wavenumbers thus 
appear to be artificially stabilised. Consequently, the selection of a maximum growth 
wavenumber determined using Saint- Venant equations should be taken as an artefact. 
This strong conclusion holds for all realistic values of (3, and H/zq. Although 

less crucial, the velocity curves also qualitatively differ in the large wavenumber limit, 
with asymptotic constant (SVSW) or growing (RANS) values of c. 
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Figure 4. Consequences of the hydrodynamical model on the dispersion relation of the linear 
instability of a fiat sediment bed. Rescaled growth rate a and phase velocity c of inclined modes 
m > as functions of ka:H, computed with either RANS or Saint- Venant stress coefficients. 
Sediment transport is assumed to be locally saturated. These curves correspond to the limit 
isat = 0. These plots have been obtained for /3/m = 31.4 (i.e. for kyH = 0.1), it,/uy-j — 2, 
= 0.1 and H/zo = 10^. The triangle indicates the location of the maximum growth rate on 
the Saint- Venant curve. The absence of such a maximum on the RANS curve shows that it is 
an artefact of SVSW approximation. The graphical convention is the same as in Fig. [2] (see 
legend) . 

The conclusion is twofold. First, the linear stability analysis previously proposed in 
the literature to explain the formation of bars suffer from the flaws of Saint- Venant 
shallow water, the observed maximum growth rate being an artefact. Second, with a 
more reliable hydrodynamical model, the growth rate does not present any maximum 
in the range of wavenumbers where bars are expected. We will now show that refining 
the sediment transport model by taking into account flux saturation transients allows 
the model to predict the formation of bars. We will see that SVSW and RANS based 
approaches can be partly reconciled when the so-called saturation length is introduced 
in the description. 

3.3. Modified dispersion relation, introducing the saturation length 

In non-homogeneous situations, the sediment flux q does not equilibrate immediately 
with the local values of the basal shear stress, but relaxes towards its saturated value 
gsat over a typical distance called the saturation length Lgat- We describe this process 
by the following first order relaxation equation: 

Lsnt (t- (f = (fsat - q, (3.9) 

This saturation equation can both describe bed load (Charru 2006, Duran et al. 2012, 
Charru et al. 2013) and suspended sediment transport (Claudin et al. 2011). In the later 
case, which is the most important for fluvial bars, the description has been validated 
using experimental data. The saturation length Lgat is an increasing function of the 
shear velocity and importantly presents a sharp increase at the transition from bed- 
load to suspended sediment transport. The saturation length is effectively set by the 
grain diameter d for bed-load (the order of magnitude is around 10 d) but rather set by 
the flow depth H for a turbulent suspension. As an example, for the Loire river shown in 
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Figure 5. Comparison of the dispersion relation of inclined modes m > at small and large 
saturation parameter Lsat/H. Growth rate as a function of kxH for (a) Laat/H = 0.1 and (b) 
Lsat/H = 10. The other parameters are P/m = 31.4 (i.e. for kyH = 0.1), Ut/u1^ — 1, F = 0.1 
and H/zo = IQ^ . The triangle (resp. circle) indicates the location of the maximum growth 
rate on the Saint- Venant (resp. RANS) curve. The most unstable modes only coincide at large 
Lsat/H. At small Lsat/H, the maximum growth rate on the SVSW curve is due to a spurious 
description of hydrodynamics at large wavenumber. The graphical convention is the same as in 
Fig. [2] (see legend). 



Fig. [T^, the saturation length increases by four orders of magnitude during floods, when 
the transition from bed load to suspended load is crossed. 

The description of the saturation transient brings a new dimensionless number in the 
problem, that we propose to call the saturation number, defined as the ratio Lsat/H be- 
tween the saturation length and the flow depth (Andreotti et al., 2012). By construction, 
the saturation number is larger than 1 for suspended load: it is typically on the order of 
10^ or 10^ during flooding events. It can also be larger than 1 in braided gravel streams 
where the particles can have a size comparable to the flow thickness. 

The saturation transient equation gives in Fourier space 

r gLf (3.10) 

The dispersion relation then becomes 

ik Q 

a - ik^c = — — — [kj; {ax + ih^) + ky {ay + iby)] . (3-11) 

1 + ikxLsat 

The new factor (1 + ifc^isat)^^ shows that the transport saturation transient has a 
stabilising effect on wavelengths comparable or smaller than Lgat ■ The dispersion relation 
is plotted in Fig. [5] for two values of Lsat/H. On all curves, the growth rate a is now 
negative at large wavenumber and presents a maximum controlled by Lgat- Comparing 
SVSW and RANS based approaches, the most unstable mode is located at different 
values of k^H for small Lg^t/H, but coincides at large L^at/H. When the saturation 
length is large enough, the transport saturation transient stabilises all wavenumbers, for 
which Saint- Venant equations are not reliable. The two predictions then coincide in the 
whole range of unstable wavenumbers. 

To conclude this paragraph, we now consider the mode m = 0, which corresponds to 
bedforms transverse to the flow. Fig. [5] compares the dispersion relation obtained for 
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Figure 6. Consequences of the hydrodynamical model on the dispersion relation of trans- 
verse modes m = 0. Rescaled growth rate g and phase velocity c as functions of kxH, com- 
puted with either RANS or Saint- Venant stress coefficients. These plots have been obtained 
for Lsat/H — 0.1, Wt/'Utjj = 1, J- = 0.1 and H/zo — 1(P . The circle indicates the location of 
the maximum growth rate on the RANS curve. In Saint- Venant approximation the instability 
of transverse modes is missed (ct < 0). The graphical convention is the same as in Fig. [2] (see 
legend) . 

SVSW and RANS hydrodynamical models, for a small but finite value of Lsat/H. As 
already noted by different authors, the Saint- Venant completely misses the instability of 
these modes, and therefore the emergence of ripples. By contrast, using a more reliable 
hydrodynamical calculation, one observes that transverse bedforms are unstable in a 
range wave- numbers going from a lower cut-off scaling on 1/ H to an upper cut-off scaling 
on 1/Lsat- We therefore reach a new important conclusion regarding the instability of a 
flat sediment bed: ripples (to ~ 0) and bars (m > 0) form by the very same instability, 
but possibly in different regimes. In order to find the most unstable mode, one must 
compare all the modes, including to = 0. Due to the artefacts of Saint- Venant approach 
at large k^H , this mode has been ignored so far in the discussion of the bar instability. 

4. Revisiting the bar instability 

In the first part of the article, we have shown that the maximum growth rate deduced 
from Saint- Venant equations, when the transport saturation transient is ignored, is a pure 
arte-fact. We have then shown that the introduction of the saturation length is necessary 
to address the bar instability, once hydrodynamics is properly described. We will now 
focus on the RANS based approach to revisit the basic properties of the instability, and 
produce stability and bedform-phase diagrams. When necessary, we will still show the 
differences with the SVSW based approach. 

4.1. Influence of the saturation parameter Lsat/H 
In order to study in more details the influence of the value of Lsat/H on the maximum 
growth rate, we have systematically varied this ratio between the two values used in 
Fig.O We display in Fig. [7] the longitudinal wavenumber fc,„ax and the wave angle ftmax 
of the most unstable mode. One observes that the saturation parameter Lsat/H triggers 
a sharp transition from transverse (to = 0) to inclined bedforms (amax switches from 
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Figure 7. Transition from transverse (m = 0) to inclined (m = 1) bedforms controlled by the 
saturation number Lsax/H. Longitudinal wavenumber fcmax (a) and angle amax (b) of the most 
unstable mode as a function of Lsat/H. The other parameters are /? = 31.4 (i.e. kyH — 0.1 for 
m = 1), Ut/u^y^ = 1, T — 0.1 and H/za = 10^. At the transition point (circles), the maximum 
growth rates of modes m = Q and ?n = 1 are equal. 



to ~ 37r/8). At the transition, located at Lsat/H ~ 1.6 for the parameters chosen to 
compute Fig. [71 the wavenumber fc^ax drops by half a decade. Below the transition, 
one observes in Fig. [7^ that k^nnxH decreases as the inverse of the saturation param- 
eter. Therefore, the wavelength of the transverse patterns (amax — 0) is proportional 
to Lsat, which corresponds to sand ripples. Beyond the transition, the selected mode 
corresponds to alternate bars (amax > 0). The most amplified wavenumber presents a 
much weaker dcpcndancc on igat- The wavelength approximately scales as iJ^/^Lg^f/C. 
Importantly, the SVSW approach completely misses this transition from transverse to 
inclined bedforms. 

4.2. Influence of the ratio u^/u^yi 
In this subsection we study the influence of the ratio u^/u^^-^ on the dispersion relation 
p. lip . We focus on the large saturation number regime, for which bars [m > 0) are more 
unstable than ripples (m — 0). This ratio enters the problem trough the flux coefficients 
Qx, bx, ay and by (see Eqs. I3.4[|3.5p . when the influence of the longitudinal and transverse 
slopes on the shear threshold is taken into account. In Figs.[8]fT0l we plot the longitudinal 
wavenumber fcmax, the wave angle amax, the growth rate (Tmax and the phase velocity Cmax 
of the most unstable mode as functions of the channel aspect ratio (3, for three different 
values of The saturation length is fixed at the value isat/-ff = 10, for which we 

have seen above that Saint- Venant and RANS based predictions are comparable. 

Close to the transport threshold, for — > 1, there is no influence of the transverse 
sediment transport (a-y — and by ^ 0). There is a stabilising effect of the longitudinal 
slope which vanishes as a tends to 7r/2. In the other limit, when Mh./u°j^ — >■ oo, the 
transverse sediment transport becomes important, but the slope effect disappears. In 
the intermediate range of shear velocities, both transverse sediment transport and slope 
effects are present. In particular, the transverse slope effect has a stabilising effect when 
a tends to tt/2. Looking only at the results based on the RANS approach, the general 
shape of the curves is very similar. For a given mode m > 0, the bar instability presents 
a threshold controlled by f3/m. In both limits of shear velocities close to the transport 
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Figure 8. Rescaled horizontal wavenumber (a), wave angle (b), growth rate (c) and phase 
velocity (d), all computed for the most unstable mode, as functions of the channel aspect ratio 
divided by the mode number. These curves have been computed for u^/uti^ — 1. The other 
parameters are Lsat/H = 10, J- — 0.1 and H/zq = 10^. The graphical convention is the same 
as in Fig. [5] (see legend) . 



threshold and far from it, this instability threshold is mostly due to a decrease of the 
stress component at large kyH (we recall that controls the destabilising effect). 
In the intermediate range of shear velocities, the transverse slope effect adds a further 
stabilising effect which leads to an increased value of the threshold aspect ratio, /3c- Above 
this threshold, one can distinguish three zones of /J/rn. At small /3/m, the wavelength 
is controlled by H/C, the angle is close to 7r/2, and the propagation speed increases 
with /3. In the intermediate range of /3/m, the angle decreases with f3; the wavelength 
roughly scales as yJHW/C, and the propagation speed still increases with /?. Finally, at 
high P /m, the growth rate drops. For a given /3, modes of higher m are therefore more 
amplified. 

Comparing the results obtained with SVSW equations to those obtained with RANS 
approach, one observes a global agreement in the trends. However, for u^/u^j^ — )■ 1 
and it*/it['jj — >■ oo. Saint- Venant equations miss the instability threshold at small /?, i.e. 




when kyH/C is too large to get correct hydrodynamical results. Furthermore, these 
equations predict a change of sign of the propagation speed Cmax, close to the transport 
threshold, which is due to a wrong sign of at large transverse wavenumbers, see 
Fig. |3^. Interestingly, the results obtained with both models are very close to each other 
in the intermediate range of shear velocities, thanks to the stabilising transverse slope 
effects: as the threshold /3c is higher. Saint- Venant equations are correct in the whole 
range of unstable modes. 

4.3. Convective-absolute instability transition 

We wish now to investigate the nature of this instability. In an absolute instability, some 
of the unstable modes propagate upstream while others propagate downstream. The 
instability then develops exponentially in time with a pattern that remains homogeneous 
in space. In a convective instability, all the unstable modes propagate in the same 
direction so that the instability develops in space and form a pattern stationary in time. 
To determine whether the bar instability is convective or absolute, the group velocity 
Vg = ^^^^^ must be analysed - sec the review of Chomaz (2005). We have computed 
Vg for the most unstable mode (cr = (imax), as well as for the marginally stable modes 




(cr = 0). There are two such marginahy stable modes and fc_, respectively above and 
below fcmax- Note that fc„ vanishes for — >• 1 (Fig. [5]). These three values of Vg are 
plotted in Fig. [TTk as a function of the channel aspect ratio. We can see that both t^max 
and w+ are positive, while v- switches from positive to negative values as (3 /m increases. 

Positive values of the group velocity v- are associated to a convective nature: all 
unstable perturbations are convected downstream while amplified. Conversely, negative 
values of v- are associated to absolute instabilities: there exist a range of k between fc_ 
and fcmax for which modes arc unstable but propagate upstream. We have plotted the 
transition between the two in Fig. [Tib in the plane /3/m vs u^/u^^ (the line for which 
V- =0). Since fc_ vanishes when u^:/u^yi ~^ 1j the convective region shrinks in this 
limit: as a consequence, an absolute instability is expected very close to this transport 
threshold. 

4.4. From alternate to multiple bars 

In paragraph 13.21 we have seen that the presence of the banks selects a discrete though 
infinite set of modes. More precisely, the boundary conditions lead to a selection of 
the transverse wavenumber: kyH = /?/m7r. Consider an initial sediment bed randomly 
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Figure 11. (a) Group velocity Vg as a function of /3/m. The three curves correspond to three 
particular points of the dispersion relation aik^H): the most unstable mode (dotted dashed line), 
and the two marginal modes at high (dashed line) and low (solid and dotted lines) wavenumbers. 
The latter changes sign, and we use a graphical convention similar to that described in Fig. O 
with a concomitant change of line symbol. These curves have been computed with Ut/uth = 2. 
(b) Stability diagram in the plane [fi /m,u* /u^]^). The line separating convective from absolute 
instabilities corresponds to vanishing group velocity. For both panels, the other parameters are 
L,at/-ff = 10, = 0.1 and H / zo = 

disturbed, say, by a white noise. Then all modes, characterised by the pair {kx,m) have 
statistically the same mean squared amplitude. The mode that emerges is the most 
unstable one, i.e. that with the maximum growth rate a with respect to kx and m. By 
contrast, in the previous paragraphs, we have analysed the properties of the maximum 
growth rate at a given transverse number m. The phase diagram of Fig. [T^] presents the 
value of m selected as a function of Lsa,t/H and /3, keeping M*/uth = 2 fixed. 

The diagram first shows that the saturation number Lgat/H controls the transition 
from transverse patterns (ripples) to bars. The second dominant feature of the diagram 
is the existence of an instability threshold. At low Lsat/-ff, the dimensionless parameter 
controlling this threshold is the ratio W/Lgat- At large Lg^t/H, the saturation length 
becomes a subdominant parameter and the bedform pattern becomes controlled by the 
channel aspect ratio /3. At small /?, a flat sand bed is stable. In an intermediate range 
of /3, alternate bars (m ~ 1) arc the most unstable modes. Finally, increasing the aspect 
ratio ^, multiple bars form with an increasing value of m. This last regime could explain 
the formation and the dynamics of bars in braided rivers (Fig. [ij)). 

4.5. Discussion and perspectives 

In conclusion, we have pointed out different artefacts in the standard theory of bar in- 
stability, which result from spurious hydrodynamical predictions of Saint- Venant shallow 
water equations at large wavenumbers. Amongst these artefacts, two must be particu- 
larly emphasised: the existence of a maximum growth rate in the linear instability and 
the stability of transverse modes. Our analysis, based instead on Reynolds averaged 
Navier- Stokes equations, reveals the fundamental importance of the relaxation of trans- 
port towards equilibrium: the dimensionless number characterising the influence of this 
dynamical mechanism is the saturation number Lgat / H turn out controls the transition 
from ripples (transverse patterns) to bars (inclined patterns). The qualitative aspects of 
the bar instability previously identified are recovered at large Lsa.t/H. The instability 
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Figure 12. Diagram showing the difTerent bedform regions in the plane {l3,Lsat /H). In the 
bottom-right dark-grey region, the bed is stable. In the left middle-grey region, ripples (m = 0) 
are the most unstable bedforms. In the remaining up-right part of the plane, higher modes 
{m > 0) are dominant: the dashed-dotted line separates alternate bars (m = 1) from double 
bars (m = 2). The dotted lines shows the transition to m = 3. The dashed line (amax ~ 7r/4) 
separates rather transverse to more elongated bars. These curves have been computed with 
u*/ut°h = 2, = 0.1 and H/zo = 10^ 

threshold and the transition from alternate bars to multiple bars are mostly controlled 
by the river aspect ratio /3, with a subdominant dependence on the saturation number. 
The instability presents at low wavenumber, a transition from upstream to downstream 
propagating modes. In particular, the so-called 'resonant conditions' where both the 
growth rate a and the propagation speed c vanish are recovered. 

Our findings have important consequences for the morphology of natural rivers. Indeed, 
small values of Lgat/H are typical of bedload transport. In this case, the primary linear 
instability of a flat erodible bed leads to the emergence of transverse ripples, i.e. to the 
mode m = 0. Conversely, large values of Lsat/H are typical of transport of particles in 
turbulent suspension (Claudin et al., 2011). Then the linear stability analysis predicts 
the emergence of alternate bars with an angle between 7r/4 and 7r/2. In conclusion, the 
transition from bed load to suspended load in sandy rivers is associated to a sharp change 
of bedforms. While centimetre-scale sand ripples form when bedload dominates, they are 
erased during floods and alternate bars then form. 

To go further in the quantitative interpretation of natural bedforms a calibration effort 
is needed, and the saturation number Lgat/^^ should be systematically measured in-situ 
and in flumes, besides standard parameters. Furthermore, as the range of unstable modes 
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is large, pattern coarsening - amongst other possible non linear effects - is expected to 
take place. As a consequence, the final wavelength measured in natural rivers or fiumes 
may be substantially different from that predicted by the linear instability analysis. 
Finally, it would be interesting to predict the evolution of the multiple bars composing 
braided streams, using the modal analysis developed here, in order to test quantitatively 
the relevance of the approach. 

We also foresee three promising theoretical perspectives on this subject. First, it 
would be interesting to derive a low dimensional model for hydrodynamics that would 
not present the artefacts of standard Saint- Venant equations (Ruyer-Quil & Manneville, 
1998, 2000; Luchini & Charru, 2010). Second, the proper asymptotic matching theory, 
including the boundary layers close to the banks, remains to be performed. Last, the 
current theory of river meandering is based on the hypothesis of weak coupling between 
bank distortions and bcdforms. with a strong emphasis on the 'resonant' conditions. It 
would be interesting to perform a complete linear stability analysis of a channel, including 
erosion and deposition in the bank region, based of RANS hydrodynamical equations. 
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